Damping and revivals of collective oscillations in a finite-temperature model of 

trapped Bose-Einstein condensation 
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We utilize a two-gas model to simulate collective os- 
cillations of a Bose-Einstein condensate at finite tempera- 
tures. The condensate is described using a generalized Gross- 
Pitaevskii equation, which is coupled to a thermal cloud 
modelled by a Monte Carlo algorithm. This allows us to 
include the collective dynamics of both the condensed and 
non-condensed components self-consistently. We simulate 
quadrupolar excitations, and measure the damping rate and 
frequency as a function of temperature. We also observe re- 
vivals in condensate oscillations at high temperatures, and 
in the thermal cloud at low temperature. Extensions of the 
model to include non-equilibrium effects and describe more 
complex phenomena are discussed. 

PACS numbers: 03.75.Fi, 05.30. Jp, 67.40. Db 



The first experimental observation of Bose-Einstein 
condensation (BEC) in magnetically trapped alkali atoms 
in 1995 p]-p| was a precursor to an explosion of interest in 
the properties of weakly-interacting Bose gases. Much of 
the subsequent theory M has focused on the dynamics of 
the condensate, including phenomena such as collective 
excitations and vortex motion. In the limit of zero tem- 
perature, one can represent the condensate by a macro- 
scopic wavefunction analogous to a classical field. In this 
case the behavior can be described in terms of the Gross- 
Pitaevskii (GP) equation, which has the form of a non- 
linear Shrodinger equation. Extension of the description 
to finite temperatures, where one must include fluctua- 
tions upon the condensate wavefunction, is a consider- 
able challenge. However, the motivation is clear, as such 
a description would allow direct comparison with exper- 
iments where a non-condensed thermal cloud is present, 
as well as revealing new phenomena such as damping of 
collective modes [J5|-[l0| and the decay of metastable vor- 
tices [0|l|]. 

Amongst the most compelling evidence for the validity 
of the GP equation at low temperatures is its quantita- 
tive agreement with experiment for low-energy collective 
modes. However, consistent theoretical descriptions at 
higher temperatures have proved far more elusive, where 
experiments have demonstrated marked frequency shifts 
and damping of the condensate modes in the presence 
of a significant non-condensed component ||[l(J . Theo- 
retical studies have tended to concentrate on one of two 
regimes, depending upon the density and temperature of 
the system. At high densities, where collisions are suffi- 
ciently rapid to force the system into local equilibrium, 



the dynamics of the condensate and thermal cloud can be 
described by a set of coupled hydrodynamical equations 
p3|"p5[ . Damping mechanisms in this case are of a dis- 
sipative type (i.e. viscosity and thermal relaxation). For 
very dilute systems or at low temperatures the mean free 
path of the elementary excitations become comparable 
to the size of the system and collisions play only a minor 
role. Damping in this collisionlcss regime is not related 
to thermalization processes but to coupling between ex- 
citations, and can be described within the framework of 
mean- field theories ||. The collisionless regime may be 
appropriate for the JILA experiments M , while the MIT 
experiments lie between the collisionless and hydrody- 
namical regimes pOj ] . 

Here we present a simple model of a finite-temperature 
BEC system, and use numerical simulations to find 
the temperature-dependent frequency and damping of 
a quadrupole mode. Essentially, we utilize a two- fluid 
approach, where the ground-state condensate and low- 
energy, highly-occupied 'classical' modes are described 
by a generalized GP equation, while the thermal cloud, 
which is composed of higher-energy excitations, is simu- 
lated using a Monte Carlo approach. One advantage of 
this model is that we do not need to invoke strong as- 
sumptions about particle collision times, so that we can 
study the intermediate region between the collisionless 
and hydrodynamical regimes. We also consistently in- 
clude the collective dynamics of the thermal cloud, which 
are particularly significant at temperatures near to the 
Bose transition. This aspect is often absent from other 
treatments, for example frequency calculations from solv- 
ing Hartree-Fock-Bogoliubov (HFB) equations Pq , p7[ . 
Finally, the model can potentially be extended to sim- 
ulate more complex situations, such as vortex decay and 
response to time-dependent probes. 

The generalized GP equation for the condensate wave- 
function, ^(r, t), is written in the Popov approximation 
(i.e. the 'anomalous' density, m — 0) as pa] 



in±*{r, t ) 



2m 



V 2 + V(r,t)+g[2n{r,t) + 



n c (r,t)]-»A(r,*)J*(r,t), (1) 

where n(r, t) is the non-condensate density, while 
n c (r,t) = |^(r, t)\ 2 is the condensate density (where the 
wavefunction is normalized to the number of condensate 
atoms, N c ). The condensate is subject to an external 



trap potential, V(r,t) = m(uj 2 r 2 + ui 2 z 2 )/2, as well as 
mean-field effects arising from the thermal cloud and 
the remainder of the condensate. These interactions are 
parameterized by the coupling constant, g = Airh 2 a/m, 
where m is the atomic mass and a is the s-wave scattering 
length. In this paper we shall concentrate on 87 Rb, where 
a = 5.5 nm. The dissipative term, A(r, t), represents col- 
lisions that transfer atoms between the condensed and 
non-condensed components. Here we shall assume local 
equilibrium between the two components, so that this 
term vanishes, A = 0. Generalization of this model to 
include this effect will be the subject of future work. 

Given the non-condensate density n as a function of 
position and time, the GP equation ([l]) can readily be 
solved using techniques discussed in our previous work 
|L8[ . Briefly, the GP equation is re-scaled for convenience, 
before being propagated over a small time-step At using 
a FFT method. The time-independent problem, corre- 
sponding to finding initial conditions of our simulations, 
can be solved by propagating in imaginary time, so that 
an arbitrary wavefunction quickly diffuses to the ground 
state solution. The equilibrium thermal distribution can 
be calculated under a semi-classical approximation p9[ , 
where in a harmonic trap of mean frequency Q the dis- 
crete energy levels are replaced by a continuous function 
£hf = P 2 /2m + V(r) + 2g[n c (r) + n(r)] — \i. This corre- 
sponds to the energy of a single particle moving within 
the mean-field. The semi-classical approximation is valid 
under the condition that ksT 3> Hu>, and when the num- 
ber of trapped atoms is large |J]. An integration over 
momentum states then simply yields 
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where z = exp[— (3{V(r)+2g(n c (r)+n(r))— fj,}] is the fu- 
gacity, At = (2Trh 2 /niksT) 1 / 2 is the thermal wavelength, 
and g a (z) = X^i z l /l a . The total number of atoms in 
the system is then given by N — N c + J dr n = N c + N, 
where N is the number of atoms in the thermal cloud. 

Self-consistent solution of (|l|) and (0) yields good 
approximations for the condensate and non-condensate 
densities at equilibrium. In particular, we iterate using 
successive evaluations of (0) and imaginary time propa- 
gation of the condensate wavefunction, to find the den- 
sities as a function of N and T. Given this initial condi- 
tion, the system dynamics resulting from a varying trap 
potential can be found by solving the time-dependent 
GP equation using propagation in real time. However, 
the problem of finding the non-condensate density be- 
comes more challenging as this is also time-dependent. 
Under the semi-classical and Hartree-Fock approxima- 
tions, and the assumption that the effective potential 
U e g(r,t) — V(r, t) + 2g \n c (r, t) + h(r, t)] varies slowly 
in space, one can show [ |l5pC| ] that the cloud of quasi- 
particle excitations may be described using a Boltzmann 
kinetic equation 
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The relationship between the phase-space distribution 
function /(r, p, i) and the non-condensate density is sim- 
ply given by 



h(r,t) 



dp 



(2Trh) : 



-f(p,r,t). 
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The right-hand term of (|3|), which provides the scatter- 
ing rate of state p, is given by an integral representing 
two-body collisions between atoms in the thermal cloud 
within the Born approximation 
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with / = /(p, r, t) and /j = /(pi,r,t). Locally, an ex- 
cited atom has the HF energy e p (r, t) — p 2 /2m+U B s(r, t). 
As above, we neglect collisions that transfer atoms be- 
tween the two components. The collision integral (|5j) 
differs from that of a classical gas [pT| by the inclusion 
of (1 + /) factors that represent Bose enhancement of 
scattering into occupied states. 

In Refs. [J13| |15|| , moments of the kinetic equation 
(0) were taken to yield hydrodynamical equations pl| . 
These can be solved explicitly under certain conditions 
using a variational method Jl5| . An alternative approach 
is to solve (||) directly. In general, this is very difficult 
owing to the six-dimensional nature of phase space. One 
possibility is to work under the assumption of sufficient 
ergodicity. Ergodicity assumes that the distribution of 
atoms in phase space depends only on their energy, e. 
Then (|2|) reduces to an equation of motion for /(e). This 
assumption is well-known in the literature and has been 
used to model evaporative cooling [^2],^3| and conden- 
sate growth J2I |26j in Bose gases. However, ergodicity 
assumes that any deformation in momentum or position 
space is isotropic, or that the ergodic mixing time is 
shorter than the elastic collision time. In general non- 
equilibrium situations this assumption is not valid. In 
addition, we are primarily interested in the gas dynam- 
ics in position space and its coupling to the condensate. 
Hence, a Monte Carlo technique B7H29] is more appro- 
priate here. In particular, we utilize a direct simulation 
Monte Carlo (DSMC) method, as performed to model 
evaporative cooling in Bose gases p8| , [29|, and described 
in detail for classical gas dynamics in ]30[|. We now dis- 
cuss our extension of this model to simulate the thermal 
cloud coupled to the condensate. 

The direct simulation Monte Carlo (DSMC) method 
was first developed by Bird to describe classical gas Sows 
p0|. It is equivalent to solving the Boltzmann equation 



in phase-space, except that it recognizes the discrete na- 
ture of the gas on a microscopic level. In principle, the 
trajectory of each atom could be followed at all times, 
so that the state of the system is completely described 
by storing (r,v) for all atoms. However, the calculation 
becomes unfeasible in the presence of interparticle colli- 
sions. Bird's method makes the key assumption that the 
free particle motion and collisions are uncoupled over a 
short time interval, At. This provides an accurate de- 
scription of the gas so long that At <C t co h, where r co n 
is the mean collision time. Hence the DSMC method 
is most appropriate for describing gases in the Knudscn 
regime, where the mean free path is much larger than the 
size of the system. The technique is therefore well suited 
to dilute alkali gases. 

First, the atoms are moved over distances appropriate 
to their velocity components, Vfe (fc G {x, y, z}), such that 
rfe+i = rfc+VfcAt, before collisions are treated. To ensure 
that collisions only take place between near neighbours, 
position space is divided into cubic cells of a size much 
smaller than the dimensions of the cloud. The number 
of atoms is counted in each cell to furnish the local den- 
sity n(r). Pairs of atoms in a cell are then chosen at 
random, and the momenta after a collision is calculated 
a priori. To account for energy and momentum conser- 
vation the collision is most conveniently treated in the 
atomic centre-of-mass frame. Two further random num- 
bers, R\ and i?2, are chosen to determine the scattering 
angles </> = 2nRi and cos# = 1 — 2R.2, where iZi 2 S [0, 1]. 
To decide whether the atoms actually do collide, the fol- 
lowing algorithm is used. First, the mean number of 
collisions locally in time At is calculated using 

p(r, t) = fi(r, t)av r At[l + /(p 3 , r, t)} [1 + /(p 4 , r, t)], (6) 

where a = 87ra 2 is the scattering cross-section for bosons 
in a hard sphere model, and v r = |v2 — Vi| is the rel- 
ative velocity of the two atoms. A 'quantum scatter- 
ing factor' is also included which represents the effect 
of Bose statistics, where P3 and P4 are the momenta of 
the collision products. To estimate the distribution func- 
tion, the number of atoms are counted within subcells 
in momentum space, which in turn are subdivisions of 
the positional cells. Strictly speaking, each phase-space 
subcell should have a volume of ft 3 , which is the mini- 
mum value allowed by the uncertainty principle. How- 
ever, the computational time required to sort the atoms 
increases linearly with the total number of subcells, and 
can become prohibitively large without some form of 
coarse graining. We therefore count the number of atoms 
A4c within larger subcells, which is renormalized to yield 
/(p, r, t) ~ A4cft 3 /V p V r , where V p and V r are the volumes 
of cells in momentum and position space respectively. We 
find that our results are largely independent of the num- 
ber of cells and subcells for sufficiently large numbers. 
For example, for the computations described below we 
use 8000 position cells subdivided into 9261 momentum 
subcells. 



As p(r,t) <§C 1, the collision probability is given by 
-Pcoii = 1 — e~^ r ' fS> . A further random number, R £ [0, 1], 
is compared to P C oii- Only if Rq < P co \\ does a collision 
takes place. Once this procedure has been repeated for 
all of the atoms in each cell, the final part of the time- 
step involves updating the atom velocities to account for 
gradients in the external potential, Vk+i = Vfc + Av^, 
where Av^ = —dkU e gAt/m. 

To simulate the coupled system, alternating Monte 
Carlo and GP propagation steps are performed during 
each time-step, At, where thermal atoms move in the 
mean- field potential, U c a. The thermal gas density is 
calculated at all points by counting atoms in each cell. 
The cells do not necessarily correspond to the GP grid, 
so cubic spline interpolation is used to smooth n. This 
prevents discontinuities in the mean-field potential, that 
may lead to instabilities in the FFT method used to prop- 
agate the condensate wavefunction. 

The first stage of the simulation is to find the ini- 
tial state for a prescribed temperature, T. The num- 
bers of condensate and thermal atoms are found by the 
semi-classical algorithm described above. The equilib- 
rium condensate density evaluated by this method, as 
well as a Maxwell-Boltzmann distribution for the thermal 
atoms, are used as an initial state for the MC-GP algo- 
rithm. The condensate is propagated through imaginary 
time while the thermal cloud relaxes to equilibrium. A 
time-dependent trap potential V(r,t) is applied and the 
system allowed to propagate in real time. 

We simulate collective excitations of N = 40000 atoms 
in a disk-shaped trap (u> r — 2n x 131 Hz, uj z — v8av). 
These are similar parameters to the JILA experiment [g| . 
We study the m = quadrupole mode, which is excited 
by a sudden 10% increase in the radial frequency, uj r 
pjj . The subsequent condensate oscillations are shown 
in Fig. ||. We clearly observe damping at higher tem- 
peratures, which is absent at T = 20 nK. Note that 
we determine the widths of both the components by 
calculating the standard deviation a\, = \J (k 2 ) — (fc) 2 , 
where we use (k n ) — Jd 3 r fc™|'I'| 2 for the condensate 
and (k n ) = J^ k™/N for the thermal cloud. To avoid 
large statistical fluctuations in the thermal component, 
especially at low temperatures, we simulate ten times the 
physical number of atoms (equivalent to repeatedly run- 
ning our simulations: a time-consuming process). For 
consistency the density and phase-space density of the 
gas are rescaled appropriately. 

We fit the condensate widths along x and y to an expo- 
nentially decaying sine function: <7k(t) = Ac~ rt coscut + 
B. The oscillation decay rate, T, and frequency, lu, are 
plotted against temperature in Fig. 0. The damping 
increases from zero at low temperatures, before tend- 
ing towards a linear dependence at intermediate values 
(T < 0.7T C °). This is in agreement with the expected be- 
haviour of Landau damping in homogeneous and trapped 
condensates, where in the limit of zero temperature the 
damping has aT~T 4 dependence, while at higher tem- 
peratures r ~ T MuMWX. We also observe quantita- 



tive agreement with previous theory [[Z| j32|j33[ ] and ex- 
periment in this regime. For example, at T = 200 nK 
{T/T° ~ 0.7) we find that V ~ 45.3 s" 1 , in fair agreement 
with the experimental value of 90 ± 40s _1 J9). Landau 
damping arises due to mean-field coupling between fluc- 
tuations in the condensate wavefunction, 5^ , and in the 
non-condensate density, 8fi [pi. Physically this is equiv- 
alent to absorption of a quantum of the collective mode 
by a thermal excitation Mm ■ We find no damping at low 
temperatures, so that Balieav damping, which is active 
even at zero temperature, is not observed. This is to 
be expected, as the mechanism involves the decay of the 
collective mode into two lower frequency excitations and 
is equivalent to coupling between 5^ and fluctuations in 
the anomalous density, Srh |8|], which are neglected in 
our model. Nevertheless, the model is still consistent be- 
cause Balieav damping is expected to be suppressed for 
the lowest modes in trapped condensates due to the dis- 
crete nature of the levels at low energy. 
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20 30 40 50 

t/ms 
FIG. 1. Quadmpole (I — 2, m = 0) oscillations of a con- 
densate at T = 20 nK (solid) and T = 160 nK (dashed). The 
width of the condensate is represented by the standard devia- 
tion along x, a x . Damping is observed at higher temperatures 
due to coupling with the non-condensed thermal cloud. 

Note that the Landau damping observed here should 
not be confused with the damping mechanism discussed 
in J14],ll5|], which is due to collisions that excite conden- 
sate atoms into the thermal cloud. For example, Landau 
damping is a three excitation process as opposed to the 
four excitations involved in collisional damping. Our ap- 
proach is justified as a first approximation because for 
relevant parameters the magnitude of collisional damp- 
ing is significantly smaller than Landau damping |34| . 

We observe a dip in the damping rate in the region 
0.7T C ° < T < 0.8T C °. This is related to an interesting 
'beating' effect in the condensate oscillation. In this tem- 
perature range the oscillations are seen to damp rapidly 
at early times, before reviving at a smaller amplitude 
after approximately ten oscillations. As a result, the fit- 
ting function tends to underestimate the damping rate. 



As shown in Fig. ||, the condensed and normal compo- 
nents oscillate at slightly different frequencies due to their 
weak coupling, and the condensate is much more highly 
damped than the thermal gas. The latter is a conse- 
quence of the more massive thermal cloud at this tem- 
perature (so that the back-action from Landau damping 
has less of an impact) and the small 'internal' damping of 
the cloud from thermalization processes. As a result the 
thermal cloud acts as a kind of energy reservoir. When 
the oscillations of the two components are in anti-phase 
the condensate oscillations are strongly damped; how- 
ever, when they are in phase the thermal cloud tends 
to drive the condensate oscillations. This beating effect 
is most noticeable when one component is significantly 
larger than the other. Correspondingly, we see the same 
effect in the thermal cloud at low temperatures. 

The condensate oscillation frequencies are also plotted 
in Fig. 0. The figure shows a small downward shift in 
frequency for T < 0.6T C ° §0,|l|,||. An increase in fre- 
quency above 0.6T° is also observed; however, this is not 
as large as that seen in the JILA experiment ||, where 
the frequency approaches 2cu x in this region. A possible 
explanation for the experimental behaviour was provided 
by Bijlsma and Stoof |35|, who suggested a cross-over be- 
tween normal modes where the two components oscillate 
in phase and anti-phase. However, as noted previously 
we find that the components oscillate at slightly different 
frequencies, and this description is inappropriate here. 
We may be able to see this effect at higher temperatures, 
though unfortunately the condensate in this regime is 
small and more sensitive to local fluctuations in the ther- 
mal cloud density, leading to unacceptable errors. The 
experimental data also suffers from large errors in this 
region, making a direct comparison difficult. 
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FIG. 2. The damping rate, F (top) and frequency w (bot- 
tom) of quadrupole m = oscillations in a cloud of 40000 
atoms. The a;-axis is plotted as function of temperature, 
T/T®, where T° = 286 nK is the ideal critical temperature, 
while the y-axes are plotted with respect to the trap frequency 
lu x = 2ii x 131 Hz. 
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FIG. 3. Quadrupolar oscillations in the thermal cloud 
(top) and condensate (bottom) at T = 240 nK. A revival of 
the condensate oscillations occurs at t ~ 40 ms. 

To summarize, we have studied frequency shifts and 
Landau damping due to mean-field coupling between the 
condensate and the thermal cloud. We also observe re- 
vivals in the condensate oscillations at high tenrperature 
due to back-coupling from the thermal cloud, illustrat- 
ing that damping is not completely irreversible. Fu- 
ture extensions to our Monte Carlo simulations could 
include collisions that excite atoms out of the conden- 
sate 1 14 15|J34|]. Along with a description of additional 
damping processes in collective excitations, this extended 
model could also be used to study condensate forma- 
tion in systems far from equilibrium. Another applica- 
tion might consider vortex dynamics at finite tempera- 
tures. A vortex in this case would be an 'obstacle' in the 
mean- field experienced by non-condensed particles, lead- 
ing to scattering and hence to a net force on the vortex. 
This should result in the expected drift of the vortex to 
the condensate edge, and allow direct determinations of 
vortex lifetimes. Similar models could also facilitate a 
fully consistent description of dissipative processes dur- 
ing probing of the condensate by a moving object p8| . 
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